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Abstract 

We compute solutions of the Lagrangian- Averaged Navier-Stokes a— model (LANS— a) for significantly 
higher Reynolds numbers (up to Re « 8300) than have previously been accomplished. This allows suf- 
ficient separation of scales to observe a Navier-Stokes inertial range followed by a second inertial range 
specific to LANS— a. Both fully helical and non-helical flows are examined, up to Reynolds numbers of 
~ 1300. The analysis of the third-order structure function scaling supports the predicted / 3 scaling; it cor- 
responds to a k~ l scaling of the energy spectrum for scales smaller than a. The energy spectrum itself 
shows a different scaling which goes as k . This latter spectrum is consistent with the absence of stretching 
in the sub-filter scales due to the Taylor frozen-in hypothesis employed as a closure in the derivation of 
LANS— a. These two scalings are conjectured to coexist in different spatial portions of the flow. The Z 3 
(E(k) ~ A;" 1 ) scaling is subdominant to k 1 in the energy spectrum, but the I 3 scaling is responsible for 
the direct energy cascade, as no cascade can result from motions with no internal degrees of freedom. We 
demonstrate verification of the prediction for the size of the LANS— a attractor resulting from this scaling. 
From this, we give a methodology either for arriving at grid-independent solutions for LANS— a, or for 
obtaining a formulation of a Large Eddy Simulation (LES) optimal in the context of the alpha models. The 
fully-converged grid-independent LANS— a may not be the best approximation to a direct numerical sim- 
ulation of the Navier-Stokes equations since the minimum error is a balance between truncation errors and 
the approximation error due to using LANS— a instead of the primitive equations. Furthermore, the small- 
scale behavior of LANS— a contributes to a reduction of flux at constant energy, leading to a shallower 
energy spectrum for large a. These small-scale features, however, do not preclude LANS— a to reproduce 
correctly the intermittency properties of the high Reynolds number flow. 

PACS numbers: 47.27.ep; 47.27.E-; 47.27.Jv; 47.50.-d 



2 



I. INTRODUCTION 



Since the degrees of freedom for high Reynolds number (Re) turbulence, as can be encountered 
in geophysical and astrophysical flows, can be very large, the implementation of their numerical 
modeling can easily exceed technological limits for computations. Furthermore, since truncation 
of the omitted scales removes important physics, e.g., of multi-scale interactions, the only ap- 
proach to a numerical study of such flows is to employ subgrid modeling of those scales. This 

nnti 

is frequently accomplished with Large Eddy Simulations (LES-see yj, |2J, yQ for recent reviews). 
This is of importance for geophysical, astrophysical and engineering applications and can have 
consequences for meteorological [Q] and climate prediction simulations [5], for instance. While 
realistic Reynolds numbers will remain out of reach for the foreseeable future, subgrid modeling 
can be an extremely useful tool in the computation of simulations for such applications. 

The incompressible Lagrangian-averaged Navier-Stokes equations (LANS— a, a— model, or 
also the viscous Camassa-Holm equation) 7, 8, 9, 10, U] is one possible subgrid model. It can 
be derived, for instance, by temporal averaging applied to Hamilton's principle (where Taylor's 
frozen-in turbulence hypothesis is applied as the closure, and also as the only approximation of the 
derivation) J12I. [13! For this reason, the momentum-conservation structure of the equations 
are retained. For scales smaller than the filter width, LANS— a reduces the steepness of steep 
gradients of the Lagrangian mean velocity and limits how thin vortex tubes become as they are 
transported (the effect on larger length scales is negligible) (SO. The a— model may also be derived 
from smoothing the transport velocity of a material loop in Kelvin's circulation theorem Jll| . Con- 
sequently, there is no attenuation of resolved circulation, which is important for many engineering 
and geophysical flows where accurate prediction of circulation is highly desirable. LANS— a 
has previously been compared to direct numerical simulations (DNS) of the Navier-Stokes equa- 
tions at modest Taylor Reynolds numbers (R x « 72 UM, R\ « 130 Q9Q, and R x « 300 Ml 611). 
LANS— a was compared to a dynamic eddy-viscosity LES in 3D isotropic turbulence under two 
different forcing functions (R\ ~ 80 and 115) and for decaying turbulence with initial conditions 
peaked at a low wavenumber (R\ « 70) and at a moderate wavenumber (R\ « 220) ]l7[|. In 
these comparisons, LANS— a was preferable in that it demonstrated correct alignment between 
eigenvectors of the subgrid stress tensor and the eigenvectors of the resolved stress tensor and 
vorticity vector. LANS— a and a related regularization, the Leray model, were contrasted with a 
dynamic mixed (similarity plus eddy-viscosity) model in a turbulent mixing shear layer (Re ~ 50) 
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19Q . LANS— a, with relatively high subfilter resolutions, was the most accurate of these three 



LES tested at this moderate Re, but it was found that the effects of numerical contamination can 
be strong enough to lose most of this potential. This could pose some limitations on its practical 
use. Quantifying those limitations is one of the goals of this present work. We will also find in this 
study that, even with sufficient subfilter resolution, LANS— a fails to represent all the neglected 
physics in a more turbulent regime (higher Re). 

The a— model also describes an incompressible second-grade non-Newtonian fluid (under a 
modified dissipation) II 1 IN - In this interpretation, a is a material parameter which measures the 
elastic response of the fluid. Either from this standpoint, from its status as a regularization of 
the Navier-Stokes equations, or, independently of any physically motivation, as a set of partial 
differential equations with proven unique regular solutions, we may analyze LANS— a without 
any LES considerations. Analyzing inertial-range scaling for LANS— a for moderate and large a, 
as well as identifying different scalings at scales larger and smaller than a is another of the goals 
of this work. In this context we also study the numerical resolution requirements to obtain well- 
resolved solutions of LANS— a (i.e., grid-independent solutions) which leads to a verification of 
the predictions of the size of the attractor in LANS— a [ill 2o|. Section HI] presents the LANS— a 
model, our numerical experiments and technique. In Section [HI] we analyze inertial-range scaling 
for LANS— a. In Section [IV] we determine the numerical resolution requirements to obtain well- 
resolved solutions of LANS— a. In Section [V] we address the LES potential of LANS— a by 
comparing a-model simulations to a 256 3 DNS (Re « 500, R x « 300), a 512 3 DNS (Re « 670, 
R x « 350), a 512 3 DNS (Re « 1300, R x « 490), a 1024 3 DNS (Re « 3300, R x « 790), and a 
2048 3 DNS (Re « 8300, R x « 1300). (The Re « 3300 simulation has been previously described 



in a study of the imprint of large-scale flows on local energy transfer H21L 12211 .) In Section [VT] 
we compare and contrast in more detail LANS— a solutions with DNS at Re ps 3300. Finally, in 
Section IVIII we summarize our results, present our conclusion, and propose future directions of 
investigation. 
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II. TECHNIQUE 



We consider the incompressible Navier-Stokes equations for a fluid with constant density, 



d t Vi + VjdjVi = -i),i> + vdjjVi + I) 

dm. = 0, 



(1) 



where v { denotes the component of the velocity field in the Xi direction, p the pressure divided 
by the density, v the kinematic viscosity, and an external force that drives the turbulence (in 
all 

k 



results, the time, t, is expressed in units of the eddy-turnover time). The LANS— a equations 
are given by 



d t Vi + UjdjVi + VjdiUj = — diTr + vdypi + F{ 

diVt = diUi = 0, 



(2) 



where Ui denotes the filtered component of the velocity field and ir the modified pressure. Filtering 
is accomplished by the application of a normalized convolution filter L : f \— > f where / is any 
scalar or vector field. By convention, we define Ui = V{. We choose as our filter the inverse of a 
Helmholtz operator, L = H^ 1 = (1 — a 2 ^) -1 . Therefore, u = g a ® v where g a is the Green's 
function for the Helmholtz operator, g a (r) = exp(— r/a) / (47ra 2 r) (i.e., the well-known Yukawa 
potential), or in Fourier space, u(k) = v(/c) /(l + a 2 k 2 ). 



We solve Eqs. CO and © using a parallel pseudospectral code Q23 



2411 in a three-dimensional 



(3D) cube with periodic boundary conditions. In most of the runs, we employ a Taylor-Green 
forcing Q, 



sin k x cos k y cos k z 
- cos kox sin koy cos koz 




(3) 



(generally, with k = 2), and employ dynamic control Il26ll to maintain a nearly constant energy 
with time. This expression Eq. © is not a solution of the Euler's equations, and as a result small 
scales are generated fast when the fluid is stirred with this forcing. The resulting flow models 
the fluid between counter-rotating cylinders [1271] and has been widely used to study turbulence, 
including studies in the context of the generation of magnetic fields through dynamo instability 
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112 811 . We also consider some runs with random and ABC Ii22l1 forcing. We define the Taylor mi- 
croscale as A = (v 2 )/{uj 2 ), and the mean velocity fluctuation as v rms = (2 J °° E(k)dk) ^ 2 . 
The Taylor microscale Reynolds number is defined by R\ = v rms \/v and the Reynolds number 
based on a unit length is Re = v rms x 1/v. 



HI. INERTIAL RANGE SCALING OF LANS— a 

A. Z 3 scaling of third-order structure function derived from the Karman-Howarth theorem for 
LANS -a 

For LANS— a, the H^(u) norm is the quadratic invariant to be identified with the energy, 



dE, 



dt 



-2uQ a , (4) 



where 



and 



Ea = ~D J D ^ U ~ a2V2u) ' Ud3x = ~~dJ d \ V ' Ud3x ' (5) 

tt a = J hjj ■ u>d 3 x. (6) 
As usual, we define the (omni-directional) spectral energy density, E a (k), from the relation 

E a = fE a (k)dadk= / E a (k)dk (7) 
Jo J Jo 

where § da represents integration over the surface of a sphere. The a— model possesses a theorem 
corresponding to the Karman-Howarth theorem j29|l for the Navier-Stokes equations and, as in the 



Navier-Stokes case, scaling of the inertial range energy spectra may be derived from it [30]. We 



summarize here the dimensional analysis argument for the LANS— a inertial range scaling that 



follows from this theorem, beginning from Equation (3.8) in Ref. 113011 . We use the short notation 
V{ = fj(x), u[ = u-(x' ,t) and r = x — x. In the statistically isotropic and homogeneous case, 
without external forces and with v = 0, taking the dot product of Eq. © with u'j we can obtain 
the equation 

d t Q l3 = -?-{T™-a 2 S™). (8) 
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The trace of this equation is the Fourier transform of the detailed energy balance for LANS— a. 



Qij = (viUj + VjUi) (9) 



is the second-order correlation tensor while 



t-r-m 

ij 



({ViUj + VjUi + VjUj + v'jU^u 171 ^ , (10) 



and 



S ij = {{d m uid i u l )u j + (drnUidjU^u'i + (g a ® T'J)vi + (g a ® r'™)^ , (1 1) 



are the third-order correlation tensors for LANS— a and r/ is the sub-filter scale stress tensor. 
For a = this reduces to the well-known relation derived by Karman and Howarth. The energy 
dissipation rate for LANS— a , e a , satisfies e a oc d t Qij. By dimensional analysis in Eq. ([8]) we 
arrive at 

1 a 2 

e a ~ -j(vu 2 + — u 3 ). (12) 



For large scales such that I 3> a, the second right hand term is ignored, u ps v, e a ~ e, and we 
arrive at the scaling of the four-fifths law, < (5v\\(l)) 3 >~ el 113 ill . Here, = [v(x + 1) — 

v(x)] • l/l is the longitudinal increment of v. The four-fifths law expresses that the third-order 
longitudinal structure function of v, S% = ((5v\\) 3 ), is given in the inertial range in terms of the 
mean energy dissipation per unit mass e by 

S v 3 = ~el, (13) 
5 

or, equivalently, that the flux of energy across scales in the inertial range is constant. We also obtain 
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34f | (hereafter, K41) energy spectrum, E[k)k ~ v 2 ~ e 2 / 3 / 2 / 3 , or, 



the Kolmogorov 1941 |3J, 
equivalently, 

E{k) ~ e 2/3 k~ 5/3 . (14) 

For small scales such that / a, however, v ~ a 2 l~ 2 u and both right hand terms are equivalent 
in Eq. (fT2l) . and our scaling law becomes 

5 3 " =< (5u||(/)) 3 >~ e a a~ 2 l 3 . (15) 



7 



Note that this scaling differs in a substantial way from the Kolmogorov scaling (~ /). For our 
small scale energy spectrum we then have 



E a {k)k ~ uv ~ a 



2/3 2/3 



(16) 



where we used u ~ a I v. The energy spectrum for scales smaller than a is then 



E a (k) ~ £^2/3^-1 



(17) 



This spectrum can also be derived from phenomenological arguments originally introduced by 
Kraichnan 113 511 . and it differs from the Navier-Stokes spectrum due to the fact that the fluid is 
advected by the smoothed velocity u which does not directly correspond to the conserved energy 
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FIG. 1: Third-order longitudinal structure function of the smoothed velocity field u, S%, versus / for large 
a LANS— q (a = 2k/3 indicated by the vertical dotted line). The scales identified with an inertial range 
are marked by vertical dashed lines and the scaling predicted by Eq. ( fT5l ). / 3 , is indicated by a solid line. 
The fitted scaling exponent Q (S^(l) ~ l^S ) is found to be Q = 2.39 ± .04. This is more consistent with 
the scaling given by Eq. (fT5T ) than K41 scaling, I 1 Eq. (fT3l) . or other proposed LANS— a scalings (indicated 
by dotted lines, see text). 



We test this prediction for LANS— a scaling at a resolution of 256 3 (y = 1.2 x 10~ 4 ) by moving 
both the forcing (k = 1) and a (k a = 2ix /a = 3) to large scales in order to increase the number 
of resolved scales for which ka > 1. In so doing, we are assuming that the scaling for large a is 
the same as for small a and large k (for evidence to this effect, see [36]). Confirmation as given 
by Eq. (TTST) is presented in Fig. Q] where we plot S% as a function of I (by convention, we plot 
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S% =< \5u\\(l)\ 3 > to reduce cancellation in the statistics). The scales identified with an inertial 
range k e [6, 10] are marked by vertical dashed lines and the predicted scaling, / 3 , is indicated by 
a solid line. We fit a scaling exponent (S%(1) ~ ) and find £3 = 2.39 ± .04. This is significantly 
steeper than the classical Kolmogorov scaling given by Eq. (fT3l ; it can thus be viewed as more 
consistent with the scaling given by Eq. ( fl3T ). It is also more consistent with I s than with other 
possible LANS— a scalings: under the assumption that the turnover time scale of eddies of size ~ / 
is determined by the unsmoothed velocity v ? we find S% (I) ~ I 5 , and if it is determined by yfy • u, 



we find £3 (/) ~ Z 4 (see, e.g., Refs. [16, 



36, 



31 



38Q). The observed scaling corresponds to none 



of these cases, and is actually closer to an evaluation of the turnover time t\ at the scale / given by 
ti ~ l/ui (with 5*3 (/) ~ I 3 ). Note that for 2D LANS-a, however, it is the case that the scaling 
is determined by the unsmoothed velocity v Q36U . We note that this is one of many differences 
between the 2D and 3D cases (e.g., ideal invariants and cascades). Another difference, which we 
shall show in Section [VH is that in 2D vorticity structures decrease in scale as a increases while 
in 3D there is a change in aspect ratio with structures getting both shorter and fatter. This may, in 
fact, be related to the shallower LANS— a energy spectrum for ka > 1 which we show in Section 
fVTl While differences are observed between the scaling shown in Fig. [Hand Eq. (fl~5l) . the error 
bars deny a K41 scaling (as well as the Z 4 and I 5 scalings) at scales smaller than a. We believe the 
discrepancy between the observed and predicted scaling can be due to lack of resolution to resolve 
properly the inertial range at sub-filter scales. We have less than a decade of inertial range and 
only 256 3 points for the statistics. As more computational resources become available, this scaling 
should be re-examined. 



B. Subdominance of the k 1 energy spectrum and rigid-body motions 

As a consequence of LANS— a's Taylor's frozen-in hypothesis closure, scales smaller than a 
can phase-lock into coherent structures and be swept along by the larger scales (see, e.g., Q30Q). 
If we assume, formally, that this "frozen-in turbulence" takes the form of "rigid bodies" in the 
smoothed velocity field (no stretching), we arrive at a much different spectrum than k~ l , Eq. 
(fTTT) . All scales smaller than a are subject to the frozen-in hypothesis and we expect to find such 
rigid bodies at these scales. We note that collections of "rigid" portions of the flow (rotating or 
non-rotating) reduce the total degrees of freedom (dof) and make physical sense with LANS— a's 
relation to second-grade fluids: these rigid bodies can be envisioned as polymerized portions of 
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FIG. 2: Spectral energy density, E{k), versus wavenumber, k, for large— a LANS— a solution. Here forcing 
(ko = 1) and a (k a = 2n/a = 3, vertical dotted line) are set at the largest scales to increase the number of 
scales for which ka > 1. Spectra are plotted for three norms: H^(u) norm (solid line), L 2 (u) norm (dotted 
line), and the L 2 (v ) norm (dashed line). As these last two norms are not quadratic invariants of LANS— a, 
we employ the H\ norm for all following results. All three spectra correspond to that derived from the 
assumption of rigid bodies in the smoothed velocity u, Eq. ( fT9l ). The vertical dashed lines are at the same 
scales as those in Fig. [TJ 

the fluid. As a matter of fact, in such structures all internal dof are frozen. These "rigid bodies" 
follow as well from the consideration of LANS— a as an initial value problem in Fourier space, 
for which we have u(k) = v(k) /(l + a 2 k 2 ). In the limit as a approaches infinity, all wavenumber 
(and spatial) dependence for v is eliminated and the entire flow is advected by a uniform velocity 
field (advection without internal degrees of freedom). 

For a rigid body there can be no stretching and, therefore, all the longitudinal velocity incre- 
ments, Sun, must be identically zero (<5u(l) = Q x 1 from basic mechanics with ft the rotation 
vector and, hence, SuM) = Su(l) ■ l/l = 0). Note that in LANS— a Eq. © the VjdiUj term 
contributes only a rotation and not a stretching of u. Such polymerization would have two con- 
sequences. Firstly, since there is no stretching, these rigid bodies would not contribute to the 
turbulent energy cascade, 

< (5mh(/)) 3 >= 0. (18) 

Secondly, the energy spectrum from dimensional analysis (u 2 ~ const, for large a /I: u = (1 + 

a 2 /l 2 )~ l v ~ l 2 v, and E a (k)k ~ uv ~ k 2 ) is 

E a (k) ~ k. (19) 
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This is, in fact, the observed LANS— a spectrum for ka 1 as is shown in Fig. |2] We verified 
that the spectrum is not the result of under-resolved runs, as is the case, e.g., in the k 2 spectrum 
observed in truncated Euler systems I139ll or in extremely under-resolved spectral simulations of 
the Navier-Stokes equations. Indeed, equipartition of the energy among all modes in a truncated 
Euler— a system should also lead to a k 2 spectrum. Along with several experiments with different 
viscosities and also with statistically homogeneous and isotropic forcing (not shown here), these 
are assurances that the observed spectrum is not a result of inadequate numerical resolution. It 
should be noted that this is the same computation for which the third-order structure function is 
shown in Fig. [Q The third-order structure function is consistent with a I 3 scaling (corresponding 
to a k" 1 energy spectrum) while the spectrum itself is k 1 . (Also shown in Fig. |2]are the L 2 {u) = 
(u 2 )/2 and the L 2 (v) = (v 2 )/2 norms which (through u ~ a 2 v/k 2 for ka ^> 1) correspond 
to k~ l and k 3 spectra, respectively. Since the analytical properties of the LANS— a solution are 
based on the energy balance, dE a /dt = —2vQ a , in the H\{u) norm, we employ this norm for all 
following results.) These two different scalings, Z 3 and k 1 , are consistent with a picture where a 
fluid has both rigid-body portions at scales smaller than a (wherein there is no turbulent cascade) 
and spatial regions between these where the cascade does take place. For the structure functions, 
a non-cascading rigid body does not contribute to the scaling and consequently the cascading 
contribution, Eq. (fT5l) . dominates. The energy spectrum, however, for the limit of k very large, is 
dominated by the k +1 term, and hence the k^ 1 component is subdominant. 

We further explore the validity of this picture by examining the spatial variation of the cubed 
longitudinal increment, (Sv\\(l)) 3 in DNS, and (5u\\(l)) 3 in LANS— a for a /I ^> 1, which in each 
case is proportional to the energy flux across a fixed scale /. (The presence of the hypothesized 
"rigid bodies" should be evident as significant portions of the flow where there is no energy flux.) 
In Fig. |3]we show visualizations of these quantities corresponding to I = 2%/ 10 (k = 10) for 
both the large-a LANS— a simulation and a highly turbulent DNS (k = 2, v = 3 x 10~ 4 ). The 
scale (k = 10) is chosen as it is in the inertial ranges of both flows. We note that for LANS— a, 
a significant portion of the flow is not contributing to the flux of energy to smaller scales (the 
filling factor for (<5m||(27t/10)) 3 < 10~ 2 is 0.67 as compared to 0.26 for the Navier-Stokes case). 
These regions can be identified as "polymerized" or "rigid bodies" in u and their locations are 
found to be robust when the I used for (5u\\(l)) 3 is varied over a factor of 2. Moreover, this is 
highlighted in the probability distribution functions (pdfs), see Fig. HI where we see the LANS— a 
pdf is more strongly concentrated around zero than the DNS. This is consistent with the idea that 
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FIG. 3: Two-dimensional slice of the cubed longitudinal increment (du\\ (27r/10)) 3 for LANS— a and 
(5v\\ (27r/10)) 3 for DNS. For all black pixels, the cubed longitudinal increment is less than 10~ 2 (approxi- 
mately consistent with rigid bodies). On the top is the large-a simulation (ko = 1, k a = 3, v = 1.2 x 10~ 4 ) 
where the filling factor (computed over the entire 3D domain) is 0.67. On the bottom is a DNS of Navier- 
Stokes (ko = 2, v = 3 x 1CP 4 ) where the filling factor is 0.26. Thus, a much greater portion of the flow is 
consistent with collections of rigid bodies for the large— a simulation. 
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FIG. 4: Pdfs of (5v\\ (2vr/10)) 3 for DNS (N = 1024, solid line), and of (Su\\ (2tt/10)) 3 for LANS— a 
(iV = 256, dashed line), and of the DNS downgraded to lower resolution (N = 256, dotted line). See Fig. 
[3]for simulation parameters. Note that both pdfs have a slight positive asymmetry consistent with a positive 
dissipation rate £r a y The LANS— a pdf is more strongly concentrated around zero consistent with the idea 
that portions of the flow (at scales smaller than a) are acting as rigid bodies. 

the internal dof of large portions of the flow (at scales smaller than a) are frozen. We point out that 
this comparison is not a LES validation, but, rather, a comparison between the dynamics of two 
different fluids at similar Reynolds numbers. One flow is a well-resolved numerical solution of the 
Navier-Stokes equations, and the other is a well-resolved solution of the LANS— a equations with 
large a. For this reason a reduced resolution (N = 256) representation for the DNS (for which 

= 1024) is not depicted in Fig. [3j We have performed such a down-sampling, however, and 
find the filling factor is reduced even more, to 0.14, and the tails of the pdf increase over the full- 
resolution analysis (dotted line in Fig. U). No inverse Helmholtz filtering, Hr l is applied to the 
DNS data. Note that this would amount to computing (6u\\ (I)) 3 in the DNS, which has no meaning 
in the dynamics of the Navier-Stokes equations (the energy flux is proportional to (8v\\ (/)) 3 ). 

We end this section with further evidence of coexistent energy spectra, k~ x and k 1 , in separate 
spatial portions of the flow. We mask out all portions of the flow that we identify with rigid 
bodies ((<5w||(27r/10)) 3 < 10" 2 , a 2D slice of which is shown in Fig. [3]). The energy spectrum 
of the remaining portion of the flow is shown in Fig. [5] as a dashed line to be compared with the 
spectrum of the entire flow shown as a solid line. The operation of spatially filtering the flow 
before computing the spectrum serves to "smear out" the energy spectrum by convolving it with 
the spectrum of the filter. Deconvolution in 3D with A" = 256 is intractable and we are, therefore, 
unable to remove this "smearing" of the energy spectrum of the cascading portions of the flow. 
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FIG. 5: Spectral energy density, E(k), versus wavenumber, k, for large— a LANS— a solution. The solid 
line indicates the spectrum as given in Fig. |2]but for a single snapshot (the same as selected for Fig. [3]). 
The dashed line indicates the spectrum wherein all portions of the flow associated with "rigid bodies" (a 
2D slice of which is shown in Fig. [3]) are removed. This provides further evidence that the flow spatially in 
between the "rigid bodies" possesses a negative power law energy spectrum (the predicted k^ 1 power law 
is shown as a solid line). 

Nonetheless, after conducting what tests we could with the filtering process (not shown here), 
we conclude that the power law of the energy spectrum of these portions is negative and, thus, 
distinctly different from that of the rigid bodies. 



IV. RESOLUTION REQUIREMENTS FOR GRID-INDEPENDENT LANS -a SOLUTIONS: 
SIZE OF ATTRACTOR 



It is useful to make a distinction between the quality of a subgrid model and effects arising from 
nonlinear interactions with discretization errors at marginal spatial resolutions (which are more 



characteristic of the discretization employed than of the subgrid model) [Jjl 4^,|41|]. Before doing 



this, we require an estimate for the total degrees of freedom for the LANS— a attractor which as 
we show, unlike for the 2D case (see [361), f° r the 3D case is reduced compared to Navier-Stokes. 
The subdominant Z 3 scaling is associated with the flux of energy to small scales and thus must 
be used to estimate the degrees of freedom of the LANS— a attractor, dof a . For dissipation the 
large wavenumbers dominate and, therefore, combining the LANS— a energy balance, Eq. ©, 
with its sub-filter scale energy spectrum, Eq. (fTTT) . allows us to implicitly specify its dissipation 
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wavenumber, k%, by 

£n ^ k 2 E a (k)dk ~ p kh^a^k-Hk ~ e a 2 / 3 « 2 / 3 (^) 2 . (20) 



Then we have, 



*S ~ <21) 



Using that the linear numerical resolution, N, must be proportional to the dissipation wavenumber 
(N > 3k") and that Re ~ v~ x , we arrive at 

AT = C k a 1/3 Re 1/2 , (22) 

or, equivalently, 

dofa = ^R e V\ (23) 
27a 



where C is an unknown constant (for further details see llllh ). We verify this prediction and 
determine the constant C through the use of a database stemming from studies in which both the 
free parameter, a (or, equivalently, k a ) and the linear resolution, N, for a set of DNS flows with 
Re « 500, 670, 1300, and 3300 are varied. In so doing, we establish the necessary numerical 
resolution for convergence to a grid-independent solution. 

Convergence to the grid-independent solution is determined by comparison of the energy spec- 
trum, E a (k), between runs with a constant filter and varying resolution. In Fig. [6£a), we make 
such a comparison for Re « 500 (N = 256 for DNS) and k a = 14 (N = 84, 96, 108, 128, and 
192 for LANS— a). We plot energy spectra compensated by k 5 ^ 3 so that a K41 k~ 5 ^ 3 spectrum 
would be flat. We see, based on comparing the energy spectra at wavenumbers smaller than k a 
to the 192 3 LANS— a spectrum, that simulations at resolutions of 96 3 and less are not converged 
while the one at 128 3 is. That is, except for the very small scales at the end of the dissipative range, 
there is very little difference between the spectra at 128 3 and at 192 3 (i.e., the solution is "grid- 
independent"). Meanwhile, for resolutions of 96 3 and less the spectra vary greatly with resolution 
(i.e., they are "unresolved"). In Fig. Ob), we collect all the results of similar studies (Re ~ 500) in 
a plot of resolution, N, versus inverse filter width, k a . (We change N for a given a, then change a 
and iterate.) Pluses correspond to grid-independent solutions, X's to under-resolved solutions, and 
squares to "undecided" runs (i.e., that are neither clearly resolved nor clearly under-resolved). The 
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FIG. 6: (Color online.) Plots for Re « 500 simulations demonstrating convergence to the grid-independent 
LANS— a solution, (a) Average energy spectra (t G [20,33], t is time in units of eddy turn-over time) 
compensated by K41 for LANS— a simulations, k a = 14: 192 3 (black solid), 84 3 (red dotted), 96 3 (green 
dashed), 108 3 (blue dash-dotted), and 128 3 (pink dash-triple-dot). The vertical dashed line denotes k a . 
Inset is a blow-up near k a where convergence can be clearly seen. LANS— a at a linear resolution of 
128 3 is approximately converged to the grid-independent solution while resolutions of 96 3 and less are 
clearly not. (b) The linear resolution of a— model simulations, N, is plotted versus k a . Simulations with 
inadequate resolution are plotted as X's, those with approximately grid-independent solutions as +'s, and 
experiments that are neither clearly resolved nor clearly unresolved as boxes. The dashed lines represent 

1/3 

N = Ckd indicating that a constant in the range 43.2 < C < 50.2 agrees with our data. This partially 
confirms the prediction of Eq. (|22l > and provides a reliable method to determine the needed resolution for a 
grid-independent LANS— a solution at a fixed Re. 

dashed lines represent Eq. (|22b with the minimal and maximal choice of C (where C = CRe 1 ^ 2 ), 
that agrees with our results (i.e., 43.2 < C < 50.2). In Fig. [7] we conduct similar studies for 
Re ~ 670. We find 49.5 < C < 51.4 and again validate the predictive power of Eq. (|22b for the 
necessary numerical resolution for grid-independent solutions. 
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FIG. 7: As Fig. Ob) but for Re ps 670 simulations. The dashed lines represent N = CkJ indicating that 
a constant in the range 49.5 < C < 51.4 agrees with our data. Note also that any power law, N oc ka, with 
0.30 < j3 < 0.46 also agrees with the data. 
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FIG. 8: Acceptable choices of C = C^Re 1 / 2 , versus Reynolds number, Re, for grid-independent LANS— a. 

1/3 

Error bars are not confidence levels, but depict the range of values consistent with our database (N = CkT) 
at the four Reynolds numbers we tested. The dashed line depicts the least-squares fit with slope 0.54 ±0.14. 
This completes the validation of Eq. (1221 ) which predicts 0.5. 

The greatest utility of the prediction, however, is due to the single constant Co which is inde- 
pendent of Reynolds number. A determination of this constant can cheaply be achieved repeating 
this process for several runs for low and moderate Re, and determines the resolution requirement 
for the highest Re attainable. The ranges of acceptable constants, C = CoRe 1 ^ 2 , for the four 
Reynolds number flows studied are plotted versus Re in Fig. [8] A power law C = Co-Re 7 fits our 
data with 7 = 0.54 ± 0.14 demonstrating the final validation of the prediction, 7 = 0.5, Eq. (|22|) . 
The value of the constant is found to be Co = 2.0 ± 0.2. We made one study for the maximally- 
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helical ABC forcing at Re « 1600 and a = 2n/2b. It is consistent with a value of C = 1.8 ± 0.1. 
We therefore conclude that the constant Co is not a strong function of the forcing employed or of 



the scale at which the system is forced. As a result, and unlike in 2D LANS— a 06H . we verify 
that the size of the attractor in 3D LANS— a is smaller than that in Navier-Stokes, which is a 
promising result if the LANS— a equation is going to be used as an LES. However, before doing 
this, an assessment of the truncation errors introduced in discretized systems (as used to solve the 
equations numerically) and a study of the optimal choice for a to capture the properties of a DNS 
is needed. We consider these problems in the following section. 

V. CAN LANS— a BE CONSIDERED AS A LARGE EDDY SIMULATION? 

In this section, we consider the LANS— a equations as a means to an end, and consider the 
solutions to their discretized equations as approximations to the Navier-Stokes solutions. We seek 
numerical approximations of LANS— a that minimize the difference to a fully resolved or direct 
numerical solution (DNS) of Navier-Stokes (i.e., we analyze the behavior of LANS— a solutions 
in the LES framework, and call here the model a "LANS— a LES", or in short "a-LES"). In the 



LES framework, LANS— en's turbulent stress tensor, r^, is given by (see, e.g., |42J]) 



r°j = H a (d k Uid k Uj + d k UidjU k - diU k djU k ). (24) 



Previous studies have not made the distinction between grid-independent LANS— a and LANS— a 



LES, though one did study convergence to grid-independent solutions at moderate Re 111911 . We 
find, however, a definite difference between the two approaches. We show in this section that, in 
fact, LANS— a combined with truncation error yields a better fit to DNS than grid-independent 
LANS— a. The resolution that yields an optimal a-LES (a terminology to be defined below) is 
also found to follow Eq. d22l) . In the Section FVAl we then address the quality and usability of the 
predictions of the LANS— a model viewed as an LES. 

A remark about nomenclature may be in order at this point. Traditionally, and for good reasons, 
LES attempt at capturing the large-scale properties of a flow with a huge Reynolds number, as 
found, e.g., in the atmosphere. In that case, the wavenumber at which the DNS is truncated 
is, at best, in the inertial range and it might even be in the energy-containing range, as for the 
atmospheric boundary layer with a Taylor Reynolds number R\ ~ 10 4 . Of a different nature are 
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the modeling methods sometimes called quasi-DNS. Here, the idea is to model a flow at a given, 
moderate Reynolds number but with an expense in computing resources lesser than if performing a 
DNS. Under-resolved DNS fall in that category; in that case, the large-scales are presumably well 
reproduced but the small scales are noisy. It is in that spirit that we now examine the properties of 
the LANS— a model. We thus qualify a model as optimal in the sense of being optimal for the class 
of LANS— a models examined herein; in order to avoid repetition, we also use the terminology of 
alpha- optimal. 



FIG. 9: (Color online.) Plot of Re « 670 simulations. Average compensated energy spectra: DNS (solid 
black line) and LANS-a simulations, k a = 41: N = 162 (red dotted), N = 192 (green dashed), and 
N = 216 (blue dash-dotted). LANS— a at a linear resolution of 192 is approximately converged to the grid- 
independent solution while a resolution of 162 is not. N = 162 does correspond, however, more closely to 
the DNS spectrum. We observe, in general, that a combination of LANS— a and truncation error yields the 
optimal a-LES. 

In Fig. [9] with k a = 41, we plot the Re ps 670 DNS spectrum (solid black line) and LANS— a 
spectra at three different resolutions. We observe that, while the N = 162 solution (dotted line, 
red online) is not converged, it is a better approximation to the DNS than the grid-independent 
LANS— a solution. For all simulations we studied, the grid-independent LANS— a solution is not 
the best approximation to the DNS. Another example is given in Fig. [101 where we plot the mean 
square spectral error normalized to make fair comparisons between large and small k a results, 



where kp is the wavenumber for the forcing scale, E{k) is the DNS spectrum (in the L 2 {v) norm), 
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E a (k) is the LANS— a spectrum (in the H^(u) norm), and n is the number of terms in the sum. 
These errors are calculated for spectra averaged over turbulent steady-state solutions: t E [16, 19] 
for Re ~ 670. We see that for a given filter or a given simulation resolution, there is a local 
minimum in the error. This minimum is a balance between truncation errors and the approximation 
error due to using LANS— a instead of the full Navier-Stokes equations. Due to these errors being, 
in some sense, in opposition, the optimal a-LES solution is found at a lower resolution than the 
grid-independent solution. Indeed, we see by examining Fig. \\0\ (a) that for a given filter the 
combination of truncation error and the LANS— a solution is a better approximation to the DNS. 
For fixed resolution, Fig. [10] (b), the optimal value for a is not zero but has some finite value. 
This local minimum error shown in the figure keeps a from going to zero (k a — > oo) in dynamical 
models 111 511 . We note, also, that the error is low for a finite range of N and k a near the minimum, 
indicating that an a-LES solution may perform well for a range of parameters near the optimal 
ones. We find the resolution for an optimal a -LES is also predicted by Eq. (122b (with C ~ 47 for 
Re ~ 670, or Co ~ 1.8). That is, optimal a-LES resolution is just below that for grid-independent 
LANS— a solutions. Having demonstrated the predictability of the resolution for grid-independent 
LANS— a and of LANS— a LES given a Reynolds number and a filter, in the following section we 
seek to determine sufficient conditions on the free parameter a for LANS— a to be a successful 
LES. 



A. Free parameter a and quality of the a-LES 

In this section, we make an analysis of the LES potential of LANS— a by considering only 
the grid-independent LANS— a solutions identified using Eq. (|22l) . Note that from the results 
discussed in the previous section, we expect LANS— a optimal grid-dependent a— LES approx- 
imations to have better performance. In the limit of a going to zero, LANS— a Eq. (0) re- 
covers the Navier-Stokes equations, Eqs. ([I]), but the question we address now is how small 
must a be for LANS— a solutions to be good approximations to Navier-Stokes solutions. There 
are several length scales that a could be related to: the forcing scale lp, the integral scale 
L = 2% E(k)k~ 1 dk/ J o °° E(k)dk, the Taylor microscale A, or the Kolmogorov dissipation 
scale t\k- Plots of the mean square spectral errors to DNS (see Eq. (125b ) versus these scales 
are shown in Fig. [TU While the general trend of errors decreasing with a is apparent in all 
cases, in Fig. [TTT a) we see a large difference between errors at varying Reynolds numbers and 
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FIG. 10: Plots for Re « 670 simulations, (a) Error (see Eq. (1251 )) versus simulation resolution for k a = 20. 
The optimal (grid-dependent) LES is for a resolution of N « 128 and has a much smaller error compared 
to the DNS than the grid-independent LANS— a solution at higher resolution, (b) Error versus k a for 
N = 128. At a given resolution the optimal value for a is not zero but occurs at a local minimal error. Any 
k a G [15, 25] has an error near the minimum indicating that an LES solution may perform well for a range 
of parameters near the optimal ones. A constant of C = CoRe 1 ^ 2 w 47 in Eq. ((22)) is found to correspond 
with optimal a-LES approximations. 

similar ratios of a to the forcing scale, If- For a linear least-squares fit, the goodness-of-fit, 
X 2 = Y.( E sq tual ~ E iq) 2 > was found t0 be X 2 = 6-2 x 10" 2 . The errors for Re w 3300 are much 
larger than for the same ratio lp/a as results at both Re ps 500 and Re ps 670. This is also the case 
for the integral scale. However, the quality of the a-LES appears to be closely tied to the ratio of 
a to the Kolmogorov dissipation scale. In Fig. [TTTb) the errors are plotted versus the ratio of the 
dissipation scale, r/ K , to a. We see a very strong dependence (x 2 — 2-5 x 10~ 2 ) between errors 
for several runs with four different Reynolds numbers indicating that the quality of the LANS— a 
LES approximation is a function of the ratio of a. to the dissipative scale. Finally, in Fig. fTTT c) 
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FIG. 11: Plot of errors, Eq. (I25I ). of grid-independent solutions compared to DNS. Asterisks are for Re « 
8300, squares for Re « 3300, triangles for Re « 670, and diamonds for Re « 500. The single right- 
most triangle in all plots corresponds to a value of a in the dissipative range (k a = 60). The norm we 
employ to measure the error, Eq. (|25T >. is no longer a good norm when dissipative scales are considered, (a) 
Errors versus lp/a. No clear correlation between LES quality and the ratio of the forcing scale to a holds 
independently of Reynolds numbers, (b) Errors versus ratio of dissipative scale, t]k, to a. The quality of 
the LES appears to be closely tied to this ratio, (c) Errors versus ratio of Taylor wavenumber, A, to a. The 
Re w 8300 experiment (asterisk) indicates that the quality of the a-LES is not tied to the Taylor scale. 
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the errors are plotted versus the ratio of the Taylor Scale, A, to a. We find x 2 = 3.1 x 1CT 2 for a 
linear least-squares fit. We note that a single experiment conducted at Re ~ 8300 (the asterisks) 
confirms that the maximal value of a is tied to the dissipation scale and not the Taylor scale. This 
is more clearly demonstrated in Fig. [12] where we plot compensated energy spectra for a nearly 
constant ratio A /a at three Reynolds numbers. We see that the maximum deviation from the DNS 
spectrum increases with Re. As A/a is the same in all cases, the optimal a is not dependent on 
the Taylor scale. 

These findings were not accessible at lower Reynolds numbers due to inadequate separation 
of scales. For example, we give in Fig. [13a) spectral flux for DNS at Re ss 500, 670, and 
3300 respectively. We define the kinetic energy transfer function, T(k), in Fourier space as 
T(k) = — J v fc ■ (<jl> x v)dV, where (■) represents the Fourier transform. For LANS— a we have 
T a (k) = — J Ufc • (u? x u)dV where u> = V x v. The flux is defined as usual from the transfer 
function as 

II (a) (fc) = / T {a) (k')dk. (26) 
Jo 

Only Re ~ 3300 (and Re ps 8300 not pictured here) demonstrates a range of nearly constant flux 
(a well-defined inertial range) before the dissipation scales. Following the scaling arguments in 
Ref. [11], one effect of the a— model is to increase the time scale for the cascade of energy to small 
scales. This reduces the flux as a increases (k a decreases) as do the hypothesized "rigid bodies;" 
this can be seen in Fig. [T3Tb). (Note that in DNS at high resolution, 80% of the flux is from local 
interactions which is strongly suppressed at scales smaller than a n2l\\ .) As dissipation dominates 
the flux for low and moderate Reynolds number, the reduced flux of the a— model has little conse- 
quence for these simulations. With a substantial inertial range, however, this reduced flux results 
in a pile-up of energy for scales larger than the dissipative scale and the spectrum approaches the 
k 1 spectrum discussed in Section [nil As a consequence of the integral conservation of energy 
(E a = J u • v) there is a corresponding decrease of energy at large scales. Consequently, as the 
inertial range increases, a must be moved to smaller and smaller scales in order for LANS— a not 
to alter scales larger than a. In summary, the a— model's reduced flux of energy to small scales is 
more crucial when the dissipation scale is farther away from a. 
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FIG. 12: Compensated averaged grid-independent energy spectra for DNS (solid) and LANS— a (dotted) 
holding the ratio of Taylor scale A to a nearly constant. Vertical dotted lines indicate k a . (a) Re « 670 
and k a = 35 (A/a = 18). (b) Re » 3300 and k a = 70 (A/a = 17). (c) Re » 8300 and fe Q = 110 
(A/a = 17). We see that the maximum deviation from the DNS increases with Re. This is due to the 
greater distance between a and the dissipative scale t}k- (Note that scales larger than k = 3 are affected by 
numerical truncation issues.) 
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FIG. 13: (Color online.) (a) Energy flux, Eq. (|26]>, for three DNS with Re « 3300 (black, solid), Re « 670 
(red, dotted), and .Re « 500 (green, dashed). No inertial range is discernible on the flux functions except 
for the highest Reynolds number case. The initial plateau followed by a bump and another plateau (for the 
case at the highest Reynolds number) is a result of the forcing employed, (b) Energy flux at Re as 3300 for 
both DNS and a— model runs; DNS is the black, solid line. See inset for LANS— a parameters. LANS— a 
gives a reduced flux which is linked to the significant pile-up of energy at high wavenumber as visible in the 
energy spectrum (see Fig [T4)) . Plots of e a versus t (not shown) also show that flux decreases (on average, at 
long times) with increasing a. 

B. Numerical savings from employing LANS— a 

If a must be directly proportional to the Kolmogorov dissipation scale, we can estimate the 
LES computational savings of the LANS— a model. For the Navier-Stokes equations we have 
dof NS oc i?e 9 / 4 and, as we verified in Section[IVj for LANS— a we have dof a = C'^k a Re 3 ^ 2 /27. 
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If k a is directly proportional to the Navier-Stokes dissipation wavenumber, k v , we arrive at 

k a « ijfc„ oc i?e 3/4 , (27) 

and, consequently, 

doi^ 5 oc i?e 9 / 4 . (28) 

Note that for free a, dof a (dof of LANS— a) is much smaller than dof^s- But, to obtain an optimal 
LES, a is tied to k v ; then the resolution requirements (dof^ ES ) are different and the decrease in 
necessary computational resolution from employing LANS— a is fixed. In fact, for the forcing and 
boundary conditions employed, we find 

dof a ES « ^dof NS . (29) 



We note that Eq. (1281) is consistent with theoretical predictions given in Ref. 



B20h. Other LES 



45Q have also exhib- 



such as the similarity model [43] and the nonlinear (or gradient) model 1144! . 
ited the characteristic that they achieve only moderate reductions in resolution and are, therefore, 
frequently used in mixed models with a Smagorinsky term (see, e.g., J3]]). That such additional 
terms will be required for LANS— a to reproduce the energy spectrum of high Re flows, may not 
be a significant factor in its usability. Note that the usual addition of extra dissipative subgrid- 
stress terms (as in the Smagorinsky model) also introduces a stronger dependence of the system 
of equations with the spatial resolution, since the filter width in such models is often associated 
to the maximum wavenumber in the box, k max . In that case, it can make more sense to use grid- 
dependent solutions of LANS— a (discussed at the beginning of Section© which give an optimal 
LANS— a LES, and can as a result give an extra gain in the computational costs. 

We also conclude that, with the scale a being tied to the dissipation scale r\ K , the model 
LANS— a behaves more like a quasi-DNS by opposition to a traditional LES. Note however that 
a factor of ~ 2.3 in resolution gain translates into a factor 27 in CPU and a factor 12 in memory 
savings, still a substantial gain. 
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FIG. 14: (Color online.) Compensated energy spectra averaged over t G [8, 9], Re « 3300. DNS is the solid 
black line and grid-independent LANS— a solutions are shown as (red online) dotted (k a = 70), (green) 
dashed (k a = 40), and (blue) dash-dotted (k a = 13) lines, respectively. A single LANS— a LES is shown 
as a (pink) dash-triple-dotted line (k a = 40, N = 384). The LES is seen to better approximate the DNS 
spectrum than the grid-independent solution for the same value of a (2-7r/40). As a is increased the energy 
spectrum approaches the k spectrum discussed in Section ITlIBi 

VI. LANS —a AT VERY HIGH REYNOLDS NUMBER 

In this section, we compare and contrast LANS— a and Navier-Stokes solutions at high 
Reynolds number. Using results of previous sections for optimal resolution and the necessary 
value of a to approximate DNS, we now evaluate both grid-independent LANS— a solutions and 
a single LANS-a LES for a highly turbulent flow (Re ps 3300, R x « 790). We calculate grid- 
independent solutions for k a = 70 (N = 512), for k a = AO (N = 512), and for k a = 13 
(N = 384). A LANS— a LES solution is computed for k a = 40 (N = 384). Averaged com- 
pensated energy spectra are shown in Fig. Q31 We see that the optimal LANS— a LES is a better 
approximation of the DNS spectra than the grid-independent LANS— a for the same value of a 
(27r/40). We also see that if a is increased further, the energy spectrum approaches the k 1 spectrum 
discussed in Section ITIIB I 

Fig [15] is a perspective volume rendering of the enstrophy density uj 2 (uj ■ uj for LANS— a) for 
the DNS, k a = 70 LANS— a, and k a = 13 LANS— a. Due to the late time depicted here (t = 9, 
longer than a Lyapunov time) there can be no point-by-point comparison between the simulations. 
However, we note that the helical structure of the vortex tubes is preserved by the a— model but 
that the tubes themselves are shorter and somewhat thicker for large values of a. As was noted 
for moderate Reynolds numbers, this is due to LANS— a suppressing vortex stretching dynamics 
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FIG. 15: (Color online.) Rendering of enstrophy density uj 2 (u> • u> for LANS— a). Due to the late time 
depicted here (t = 9, longer than a Lyapunov time) there can be no point-by-point comparison between 
the simulations. Instead, regions with approximately the same dimensions are selected around vortex tubes. 
Velocity v field lines are also shown illustrating the helical nature of the tubes which is seen to be captured 
by LANS— a. (a) DNS. The thick bars represent, from top to bottom, the Taylor scale A and the dissipative 
scale rjK, respectively. For LANS— a results the scale a is depicted between these two. (b) k a = 70, 
N = 512. (c) k a = 13, N = 384. We see that, for large values of a, the vortex tubes become shorter and 
somewhat thicker. 

without changing its qualitative features [|90. This is in contrast to 2D LANS— a where the vorticity 
structures are seen to get thinner as a increases n36fl . This could also be related to the scaling 
differences between 2D and 3D LANS— a. It has been claimed that the development of helical 
structures in turbulent flows can lead to the depletion of nonlinearity and the quenching of local 
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interactions [4 



3 



4711 . The depletion of energy transfer due to local interactions at some cutoff in 
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Consistent with 



wavenumber is also believed to bring about the bottleneck effect II22L 
these results, in 2D LANS— a (where the vorticity structures are more fine than Navier-Stokes) the 
spectrum is steeper and in 3D LANS— a (where the vorticity structures are shorter but fatter than 
Navier-Stokes) the spectrum is shallower. 




FIG. 16: (Color online.) Compensated 3rd-order structure function versus length I (a horizontal line scales 
with I). Structure functions corresponding to the Karman-Howarth theorem are depicted (©3 for DNS, 
©3 = ((5u) 2 5v) for LANS— a). Labels are as in Fig. [14] The dotted vertical lines indicate the various a's. 
A small inertial range for the DNS near I = 1 is reproduced by LANS— a. The largest a (2-7r/13) exhibits 
a second inertial range at scales just smaller than a (((Su) 2 8v) ~ I is consistent with Eq. (fT5T>). 

Figure [16] shows the third-order (mixed) structure functions corresponding to the Karman- 
Howarth theorems versus length /. For the DNS, we show 63 = (5v 3 ) and 63 = ((5u) 2 Sv) 
for LANS— a. The dotted vertical lines indicate the various a's. A small inertial range for the 
DNS near / = 1 is reproduced by all LANS— a results. The largest a (2%/13) exhibits a second 
inertial range at scales just smaller than a (((5u) 2 5v) ~ I is consistent with Eq. (fT"5T)). We note 
this is the first demonstration of third-order structure functions in LANS— a consistent with a K41 
inertial range followed by an a inertial range and finally a dissipative range. Next, we observe the 
scaling of the longitudinal structure functions, 

S p (l) = (\6v l{ \ p ), (30) 



29 



where we again replace the if * norm for the L 2 norm in the case of LANS— a, 



S"(l) = (\SuM p / 2 ). 



(31) 



We utilize the extended self- similarity (ESS) hypothesis [5l|, [52J, |53|] which proposes the scaling 



S p (l) oc S 3 (l)t» 



(32) 



or, for LANS -a, 

S%(t) oc ((6u) 2 6v)*». (33) 

We display our results in Fig. [T7j We note that for LANS— a, the third-order exponent is not equal 
to unity, contrary to the Navier-Stokes case. The Karman-Howarth theorem implies ((5u) 2 5v) ~ /, 
not (I) ~ /. We measured the deviation from linearity for each experiment (not depicted here) 
and found that LANS— a becomes more intermittent as a increases (k a = 13 is slightly more 
intermittent than the DNS). As artificially dropping local small-scale interactions gives enhanced 
intermittency J54, 55], this increased intermittency is the expected result of LANS— a reducing in- 
teractions at scales smaller than a. We note, however, that even with such a large filter, LANS— a 
is a good approximation to the intermittency properties of the DNS. This is surprising given its en- 
ergy spectrum and reduced flux in the inertial range. It is probably linked to the fact that LANS— a 
preserves global properties (in an H 1 sense) of the Navier-Stokes equations and that these proper- 
ties are important to the dynamics of small scales as measured by high-order structure functions. 



VII. CONCLUSIONS 



We computed solutions of the Lagrangian- Averaged Navier-Stokes a— model (LANS— a) in 
three dimensions for significantly higher Reynolds numbers (up to Re ~ 8300) than have previ- 
ously been accomplished and performed numerous forced turbulence simulations of LANS— a to 
study their equilibrium states. The results were compared to DNS for Re ~ 500, 670, 3300, and 
8300 , the last performed on a grid of 2048 3 points. We note that there are two ways to view the 
LANS— a simulations: as converged or "grid-independent" solutions of the LANS— a equations 
or as large-eddy simulations (a— LES) which include grid effects. We found a definite difference 
between the two approaches in that the fully-converged grid-independent LANS— a is not the best 
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FIG. 17: (Color online.) Structure function scaling exponent £ p versus order p. Black X's are shown 
for the DNS. Grid-independent LANS— a are shown as (red online) boxes (k a = 70), as (green) triangles 
(k a = 40), as (blue) diamonds (k a = 13). LANS— a LES (k a = 40, N = 384) is shown as (pink) asterisks. 
The dashed line indicates K41 scaling and the solid line the She-Leveque (SL) formula 1156(1 . 



approximation to a DNS of Navier-Stokes. Instead, the minimum error is a balance between trun- 
cation errors and the approximation error due to using LANS— a instead of the full Navier-Stokes 
equations. Due to these errors being, in some sense, in opposition, the optimal a-LES solution 
was found at a lower resolution than the grid-independent solution (the error was low for a finite 
range of N and a near the minimum, indicating that a LANS— a viewed as an LES solution may 
perform well for a range of parameters). Unlike the 2D case Il36ll . 3D LANS— a has been shown 
to be a subgrid model (i.e., it reduces the resolution requirements of a given computation). This 
difference between 2D and 3D LANS— a indicates that other a— models (as the LAMHD— a Eqs. 



a 



5811 or the BV— a Eqs. Il42h ) may behave differently and studies of these systems at high 



resolution may be required. 

We confirm the presence of the theoretically predicted I 3 sc aling of the third-order structure 
function (corresponding to a A; -1 scaling of the energy spectrum) lllll. ll6l. l37n through its bound on 



the number of degrees of freedom for LANS— a lllih . in the structure functions of the smoothed 
velocity in simulations with large a, and in the spectrum of specific spatial portions of the flow. 
In so doing, we have validated the predictive power of the bound dof a < Ca~ l Re 3 ^ 2 , for the nu- 
merical resolution for grid-independent LANS— a solutions and for optimal LANS— a LES (with 
a separate constant of proportionality). The great utility of the prediction is that the single con- 
stant can cheaply be determined at low and moderate Reynolds number and predicts the resolution 
requirement for the highest Reynolds numbers attainable. We further found no great change in 
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this single constant when employing the non-helical Taylor-Green or the maximally-helical ABC 
forcings. 

However, the small scale (ka ^> 1) LANS— a spectrum was observed to be k +1 . We attribute 
this to the frozen-in-turbulence closure employed in deriving the a— model. For scales smaller 
than a, portions of the smoothed flow u are locked into "rigid bodies." By "rigid bodies," we 
mean the internal degrees of freedom are frozen and these portions give no contribution to the 
energy cascade. This is consistent both with the observed k +1 spectrum and with field increments 
5u\\ being observed to be approximately zero over a large portion (compared to Navier-Stokes) 
of the flow. The turbulent energy cascade occurs in the space between these "rigid" portions. 
While the A;" 1 portions are subdominant to the k +l portions in the energy spectrum, they prevail 
in the cascade and hence both the structure functions and the degrees of freedom of the LANS— a 
attractor. 

We find that both of these scalings (k +1 and k~ l ) contribute to a reduction of flux at constant 
energy (i.e., the dissipation is reduced as has previously been observed in 2D calculations |p9|D . 
This leads to a shallower (or even growing) energy spectrum as a increases. Thus, for LANS— a 
viewed as an LES to reproduce the Navier-Stokes energy spectrum it is necessary that a be not 
much larger than the dissipation scale (a ;$ Ar] K independent of Reynolds number); in that sense, 
it can be considered as a quasi-DNS as opposed to a traditional LES, substantially larger Reynolds 
numbers being modeled in the latter case, leading to substantially larger gain in resolution. As 
a consequence, the computational savings of LANS— a is fixed and not a function of Reynolds 
number. (However, and unlike the 2D case, the 3D a— model does give a computational saving 
when used as a LES.) This result was not accessible at lower Reynolds numbers due to inadequate 
separation of scales. However, in one previous study for decaying turbulence with energy initially 
mostly at low wavenumbers (k = 3), it was evident that as time evolved and energy moved to 
smaller scales, the resolution requirements of LANS— a increased II 1711 . Other LES such as the 
similarity model [43!] and the nonlinear (or gradient) model [44, 45I have also exhibited the char- 
acteristic that resolution may be decreased only modestly and are, therefore, frequently used in 
mixed models with a Smagorinsky term (see e.g., [|3|]). That such additional terms will be required 
for LANS— a to reproduce the energy spectrum of high Re flows, may not be a significant factor 
in its usability. 

We compared and contrasted LANS— a to a DNS at Re ~ 3300 considering both structures 
and high-order statistics such as the longitudinal structure functions which are related with inter- 
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mittency. With an appropriate choice of a we were able to observe a Navier-Stokes inertial range 
followed by LANS— a inertial range at scales smaller than a. For this second inertial range we 
again observed a k +1 energy spectrum. As a increased, we noted a change in the aspect ratio of 
vortex tubes (they became shorter and fatter). This can be related to quenching of local small- 

sca 



(22L 



e interactions at scales smaller than a and, thus, to the shallower spectrum for 3D LANS— a 
Therefore, in 2D LANS— a (where the vorticity structures are more fine 



than Navier-Stokes) the spectrum is steeper []36U and in 3D LANS— a (where the vorticity struc- 
tures are shorter but fatter than Navier-Stokes) the spectrum is shallower. Finally, an examination 
of the longitudinal structure functions indicate that intermittency is increased as the parameter a 
is increased consistent with the suppression of local small-scale interactions at scales smaller than 



a [54, 



5J. 



The elimination of the faster and faster interactions among smaller and smaller scales through 
the modified nonlinearity in LANS— a (together with the discrepancy between its solutions and 
Navier-Stokes solutions) highlights the importance of these interactions down to scales only 
slightly larger than the dissipative scale. That is, by removing these interactions anywhere in 
the inertial range (e.g., a <; 4^), the resulting energy spectrum was found to differ from the DNS 
at scales larger than a. The intermittency properties of the DNS, however, were well reproduced 
even with large filters. Noting this, if LANS— a 's k 1 energy spectrum is not important for a given 
application, much greater reductions in resolution can be achieved. Future work should address 
whether this may be remedied in a LANS— a LES by the inclusion of another (dissipative) model 



58QW 



lether this problem is 



61 



62D. The effect of 



for these interactions, or (in the case of magneto-hydrodynamics [|57L 
less significant because of the presence of greater spectral nonlocality [|60l . 
LANS— a on the detailed scale-by- scale energy transfer should also be investigated as our results 
indicate that a model for local small-scale interactions would improve the a— model. Another 
direction of future research is to explore other reduced LANS— a models, Clark— a and Leray-a, 
which break the frozen-in-turbulence closure and, also, the conservation of circulation. Finally, 
note that because of its greater mathematical tractability, LANS— a possibly allows for a better un- 
derstanding of multi-scale interactions in turbulent flows thus modeled; therefore, detailed studies 
such as the one presented here may, in fine, allow for a better understanding of turbulence itself. 
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